home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / complex.c < prev    next >
C/C++ Source or Header  |  1990-10-02  |  4KB  |  246 lines
  1. /* complex - Complex number functions                                  */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlproto.h"
  11. #include "xlsproto.h"
  12. #include "osproto.h"
  13. #else
  14. #include "xlfun.h"
  15. #include "xlsfun.h"
  16. #include "osfun.h"
  17. #endif ANSI
  18. #include "xlvar.h"
  19.  
  20. Complex makecomplex(x)
  21.     LVAL x;
  22. {
  23.   Complex c;
  24.   
  25.   if (numberp(x)) {
  26.     c.real = makedouble(x);
  27.     c.imag = 0.0;
  28.   }
  29.   else if (complexp(x)) {
  30.     c.real = makedouble(realpart(x));
  31.     c.imag = makedouble(imagpart(x));
  32.   }
  33.   else xlerror("not a number", x);
  34.   return(c);
  35. }
  36.  
  37. double phase(c)
  38.     Complex c;
  39. {
  40.   double phi;
  41.   
  42.   if (c.real == 0.0) {
  43.     if (c.imag > 0.0) phi = PI / 2;
  44.     else if (c.imag == 0.0) phi = 0.0;
  45.     else phi = -PI / 2;
  46.   }
  47.   else {
  48.     phi = atan(c.imag / c.real);
  49.     if (c.real < 0.0) {
  50.       if (c.imag > 0.0) phi += PI;
  51.       else if (c.imag < 0.0) phi -= PI;
  52.       else phi = PI;
  53.     }
  54.   }
  55.   return(phi);
  56. }
  57.  
  58. double modulus(c)
  59.     Complex c;
  60. {
  61.   return(sqrt(c.real * c.real + c.imag * c.imag));
  62. }
  63.  
  64. Complex cart2complex(real, imag)
  65.     double real, imag;
  66. {
  67.   Complex val;
  68.   val.real = real;
  69.   val.imag = imag;
  70.   return(val);
  71. }
  72.  
  73. LVAL cvcomplex(c)
  74.     Complex c;
  75. {
  76.   return(newdcomplex(c.real, c.imag));
  77. }
  78.  
  79. Complex polar2complex(mod, phi)
  80.     double mod, phi;
  81. {
  82.   Complex val;
  83.   double cs, sn;
  84.   
  85.   if (phi == 0) {
  86.     cs = 1.0;
  87.     sn = 0.0;
  88.   }
  89.   else if (phi == PI / 2) {
  90.     cs = 0.0;
  91.     sn = 1.0;
  92.   }
  93.   else if (phi == PI) {
  94.     cs = -1.0;
  95.     sn = 0.0;
  96.   }
  97.   else if (phi == -PI / 2) {
  98.     cs = 0.0;
  99.     sn = -1.0;
  100.   }
  101.   else {
  102.     cs = cos(phi);
  103.     sn = sin(phi);
  104.   }
  105.   val.real = mod * cs;
  106.   val.imag = mod * sn;
  107.   return(val);
  108. }
  109.  
  110. Complex csqrt(c)
  111.     Complex c;
  112. {
  113.   return(polar2complex(sqrt(modulus(c)), phase(c) / 2));
  114. }
  115.  
  116. Complex cexp(c)
  117.     Complex c;
  118. {
  119.   return(polar2complex(f_exp(c.real), c.imag));
  120. }
  121.  
  122. Complex clog(c)
  123.     Complex c;
  124. {
  125.   double mod;
  126.   
  127.   mod = modulus(c);
  128.   checkfzero(mod);
  129.   return(cart2complex(f_log(mod), phase(c)));
  130. }
  131.  
  132. Complex cexpt(cb, cp)
  133.     Complex cb, cp;
  134. {
  135.   if (modulus(cp) == 0.0) return(cart2complex(1.0, 0.0));
  136.   else  return(cexp(cmul(clog(cb), cp)));
  137. }
  138.  
  139. Complex cadd(c1, c2)
  140.     Complex c1, c2;
  141. {
  142.   return(cart2complex(c1.real + c2.real, c1.imag + c2.imag));
  143. }
  144.  
  145. Complex csub(c1, c2)
  146.     Complex c1, c2;
  147. {
  148.   return(cart2complex(c1.real - c2.real, c1.imag - c2.imag));
  149. }
  150.  
  151. Complex cmul(c1, c2)
  152.     Complex c1, c2;
  153. {
  154.   double m1, m2, p1, p2;
  155.   
  156.   m1 = modulus(c1);
  157.   p1 = phase(c1);
  158.   m2 = modulus(c2);
  159.   p2 = phase(c2);
  160.   return(polar2complex(m1 * m2, p1 + p2));
  161. }
  162.  
  163. Complex cdiv(c1, c2)
  164.     Complex c1, c2;
  165. {
  166.   double m1, m2, p1, p2;
  167.   
  168.   m1 = modulus(c1);
  169.   p1 = phase(c1);
  170.   m2 = modulus(c2);
  171.   p2 = phase(c2);
  172.   checkfzero(m2);
  173.   return(polar2complex(m1 / m2, p1 - p2));
  174. }
  175.  
  176. Complex csin(c)
  177.     Complex c;
  178. {
  179.   Complex x1, x2, val;
  180.   
  181.   x1 = cart2complex(-c.imag, c.real);
  182.   x2 = cart2complex(c.imag, -c.real);
  183.   val = csub(cexp(x1), cexp(x2));
  184.   return(cart2complex(val.imag / 2.0, -val.real / 2.0));
  185. }
  186.  
  187. Complex ccos(c)
  188.     Complex c;
  189. {
  190.   Complex x1, x2, val;
  191.   
  192.   x1 = cart2complex(-c.imag, c.real);
  193.   x2 = cart2complex(c.imag, -c.real);
  194.   val = cadd(cexp(x1), cexp(x2));
  195.   return(cart2complex(val.real / 2.0, val.imag / 2.0));
  196. }
  197.  
  198. Complex ctan(c)
  199.     Complex c;
  200. {
  201.   Complex e1, e2, val;
  202.   
  203.   e1 = cexp(cart2complex(-c.imag, c.real));
  204.   e2 = cexp(cart2complex(c.imag, -c.real));
  205.   val = cdiv(csub(e1, e2), cadd(e1, e2));
  206.   return(cart2complex(val.imag, -val.real));
  207. }
  208.  
  209. Complex casin(c)
  210.     Complex c;
  211. {
  212.   Complex sx, ix, val;
  213.   
  214.   sx = cmul(c, c);
  215.   sx = csqrt(cart2complex(1.0 - sx.real, - sx.imag));
  216.   ix = cart2complex(-c.imag, c.real);
  217.   val = clog(cadd(ix, sx));
  218.   return(cart2complex(val.imag, -val.real));
  219. }
  220.  
  221. Complex cacos(c)
  222.     Complex c;
  223. {
  224.   Complex sx, val;
  225.   
  226.   sx = cmul(c, c);
  227.   sx = csqrt(cart2complex(1.0 - sx.real, - sx.imag));
  228.   sx = cart2complex(-sx.imag, sx.real);
  229.   val = clog(cadd(c, sx));
  230.   return(cart2complex(val.imag, -val.real));
  231. }
  232.  
  233. Complex catan(c)
  234.     Complex c;
  235. {
  236.   Complex sx, ix, val, one;
  237.   
  238.   sx = cmul(c, c);
  239.   sx = cart2complex(1.0 + sx.real, sx.imag);
  240.   one = cart2complex(1.0, 0.0);
  241.   sx = csqrt(cdiv(one, sx));
  242.   ix = cadd(one, cart2complex(-c.imag, c.real));
  243.   val = clog(cmul(ix, sx));
  244.   return(cart2complex(val.imag, -val.real));
  245. }
  246.